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We predict a new mechanism of surface plasmon amplification in graphene-insulator-graphene 
van der Waals heterostructures. The amplification occurs upon the stimulated interlayer electron 
tunneling accompanied by the emission of a coherent plasmon. The quantum-mechanical calculations 
of the non-local high-frequency tunnel conductivity show that a relative smallness of the tunneling 
exponent can be compensated by a strong resonance due to the enhanced tunneling between electron 
states with collinear momenta in the neighboring graphene layers. With the optimal selection of the 
barrier layer, the surface plasmon gain due to the inelastic tunneling can compensate or even exceed 
the loss due to both Drude and interband absorption. The tunneling emission of the surface plasmons 
is robust against a slight twist of the graphene layers and might explain the electroluminescence 
from the tunnel-coupled graphene layers observed in the recent experiments. 


The ultrarelativistic nature of electrons in graphene 
gives rise to the uncommon properties of their collective 
excitations - surface plasmons [iHsI . The deep subwave¬ 
length confinement (3|, the unconventional density de¬ 
pendence of frequency 0,i , and the absence of Landau 
damping 0] are probably the most well-known features 
of plasmons in graphene-based heterostructures. Among 
more sophisticated predictions there stand the existence 
of weakly damped transverse electric plasmons 0| and 
quasi-neutral electron-hole sound waves near the neutral¬ 
ity point 0, Si- Some peculiar types of plasmons can 
be excited in the graphene p — n junctions 0|, field- ef¬ 
fect transistors [l3, El , optoelectronic modulators [l^ , 
and nanomechanical resonators [l^ engaging for the im¬ 
proved device performance at the terahertz frequencies. 


Unfortunately, the experimental studies of graphene 
plasmons are yet unable to confirm or refute many of 
these predictions. To achieve an extreme plasmon con¬ 
finement, one has to sacrifice their propagation length. 
The latter is of the order of several micrometers at the 
infrared frequencies and is limited by the inter¬ 

band absorption in intrinsic samples and somewhat lower 
Drude absorption in the doped ones [l^. The experi¬ 
mentally reported quality factors of grap hene plasmons 
reach only five for graphene on Si02 and 25 for 

graphene encapsulated in hexagonal boron nitride 17| at 
room temperature. In the latter case, the damping is due 
to the scattering by the intrinsic acoustic phonons [l^ 
and can be suppressed only by lowering the temperature. 


Instead of reducing the plasmon loss, it is possi¬ 
ble to overcome the damping by introducing the gain 
medium which can replenish the energy being dissi¬ 
pated upon scattering. This idea has stimulated the 
re-examination of various ’classical’ plasma instabilities 
in graphene, including the beam and resistive instabili¬ 


ties [ 3 , Dyakonov-Shur self-excitation 13 , 111 : and gen¬ 
eration due to the negative differential conductance jl9j . 
On the other hand, the plasmon gain can be provided 
by the photogenerated electrons and holes recombining 
with plasmon emission , which opens up the prospects 
of graphene-based spasers [2l|. However, relatively fast 
nonradiative recombination in graphene [22l . l23| hinders 
the practical implementation of those structures. 

In this paper, we study the inelastic tunneling in 
graphene-dielectric-graphene van der Waals heterostruc¬ 
tures accompanied by the emission of surface plasmons 
(SPs), and the possibility of the coherent SP amplifica¬ 
tion due to the tunneling gain. The inelastic plasmon- 
assisted tunneling was first observed in the late 70’s; it 
was shown to be responsible for the light emission from 
metal-insulator-metal tunneling diodes 2^. Afterwards, 
the tunneling excitation of SPs was demonstrated using 
the scanning probes above the metal surfaces 25-2^. 
However, the possibility of the SP gain due to the tunnel¬ 
ing to exceed the loss has been never considered realistic 
due to the smallness of the tunnel exponent, thus all ex¬ 
periments to date have reported only on the spontaneous 
SP emission. To achieve large tunneling gain, one needs 
some resonant feature to compensate the smallness of the 
tunnel exponent. As an example, such a resonance occurs 
in the quantum cascade lasers due to the alignment of 
energy subbands in the neighboring quantum wells (23 - 
[snj l. In this paper, we show that the frequency depen¬ 
dence of plasmon tunneling gain in double graphene layer 
structures also exhibits a strong resonance due to the en¬ 
hanced interaction between the Dirac electrons having 
collinear momenta 3]|. In twisted graphene layers, the 
finite momentum of emitted plasmon can bridge the elec¬ 
tron dispersions in the neighboring layers together [^ . 
which makes the effect of plasmon tunneling emission ro- 
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bust against a slight layer twisting. We also discuss the 
recent experimental observations of the terahertz emis¬ 
sion from the double graphene layer heterostructures [33 1 
and address the role of the inelastic-tunneling plasmon 
excitation in the observed spectra. 





FIG. 1. (a) Schematic views of the double-graphene-layer tun¬ 
nel heterostructure, (b) The spectra of the supported plasmon 
modes calculated for d = 3nm, ef = 100 meV, T = 300 K. 
Insets in show the spatial profiles of the plasmon potential 
in the optical (top) and acoustic (bottom) modes, (c) The 
band diagram of the double graphene layer hetero-structure. 
The straight red arrows correspond to electron tunneling be¬ 
tween graphene layers accompanied by the emission of surface 
plasmons (wavy arrows). 


The double-graphene layer heterostructures shown 
schematically in Fig. 1 a support optical and acous¬ 
tic surface plasmon modes having symmetric and anti¬ 
symmetric distributions of electric potential [sj, 1^ . 
Their spectra are studied in detail in the absence of tun¬ 
neling: the dispersion of symmetric (optical) mode is 
square-root, uj+ ~ where vg is the Fermi 

velocity in graphene, qp is the Fermi wave vector, ac = 
e^/uhvo is the coupling constant, while the dispersion of 
the antisymmetric (acoustic) mode is sound-like, = sq 
(Fig. lb). Its velocity approaches the Fermi velocity vq 
as the interlayer distance d decreases, but newer falls be¬ 
low it into the region of Landau damping. In realistic 
tunnel-coupled double layers {d fv 3 nm, Sp ~ 0.1 eV), 
s = I.ISuq. Just the acoustic mode has a nonzero aver¬ 
age field strength between the two layers and thus can 
induce interlayer tunneling [^ . a process in which an 
electron passes from one layer to another with emission 
of a coherent plasmon (Fig. 1 c). The tunneling can be 
included into the dispersion law of acoustic SPs by con¬ 
sidering the tunnel current density Jpquj as a generation- 
recombination term in the continuity equations for elec¬ 
trons in a single layer 37|. In the linear response mode, 
this current is proportional to the interlayer potential 
difference, - ‘Pb.^q), where Gq^ is the 

tunnel conductance. In these notations, the dispersion 
law of the acoustic SP mode becomes 


iujK f qd\ 2Gq^ 

— (l + COth-j+,T,„ + — = 0, 


( 1 ) 


where aqui is the in-plane conductivity of graphene. Its 
real part, RetTqi.,;, is positive and leads to the SP damp¬ 
ing. At the same time, ReGqtj can be negative, at least, 
at low energy of SP quanta hw. The negative tunnel 
conductance appears as the number of electrons with the 
energy E in top layer exceeds the number of electrons 
with energy E — fyjj hi the bottom layer; such carrier 
distribution can be viewed as a ’remote population in¬ 
version’. The sign of the quantity Re \<7q^ + 2Gqtj/g^] 
thus determines whether the SPs propagating along the 
double-layer are being damped or amplified. 

An only missing ingredient to judge on the possibil¬ 
ity of the net plasmon gain is the theory of non-local 
[q ^ 0) high-frequency tunnel conductivity of van der 
Waals structures; the latter, to the best of our knowl¬ 
edge, has been studied only in the DC limit [H, 1^ . 
To construct this theory, we start with the tight-binding 
Hamiltonian of the tunnel coupled layers in the presence 
of the propagating plasmon (we set = 1): 


H = 


(Hg+ 

V f* 


f 

Hg- 


+ Llint = Hg + Lfint. 


( 2 ) 


Here, Hg± = vo<^P ± /A/2 describe separate graphene 
layers, A is the energy spacing between the Dirac points, 
p is the in-plane momentum operator, / is the iden¬ 
tity matrix, and E is the tunneling matrix. For the 
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sake of analytical traceability, we choose T in its sim¬ 
plest form which is applicable for AA-stacking of aligned 
graphene layers [s^, , T = 07, where 72 is the ’tun¬ 

neling frequency’. The interaction part, Hinti describes 
the electron-plasmon coupling; its matrix elements are 
calculated as the overlap of the eigen functions of Hq 
with the electron potential energy in the field of SP 
eip{x,z,t) = where (^o is the amplitude 

of potential on the top layer and s{z) is the dimensionless 
’shape function’ (see Fig. 1 b). 

A relatively simple form of the 7^-matrix allows us to 
treat the tunneling non-perturbatively The good 

quantum numbers of the eigen states of Hq are the in¬ 
plane electron momentum p, the band index s = ±1 (s = 
-|-1 for the conduction and s = — 1 for the valence band), 
and an extra index / = ±1 governing the z-localization of 
electron wave function. At large bias voltage, I = -1-1 and 
I = —1 correspond to the states localized on the top and 
bottom layers, respectively, while at small bias I = -1-1 
and I = —1 describe the anti-symmetric and symmetric 
states of the tunnel-coupled quantum wells. The eigen 
energies of these states are 


e‘^ = svop + Wn‘^ + AyA. 


( 3 ) 


The evaluation of tunnel current density is based on 


the relation Jxqw 


^Tr 


J±c 




p+q,p' 


, where 


is the electron density matrix calculated up to the first 
order in electron-plasmon interaction, g = 4 is the spin- 
valley degeneracy factor, and the indices a = {p, 2,s}, 
/3 = {p',F,s'} run over all quantum numbers. The ex¬ 
plicit form of the tunnel current operator is found from 
the ’equation of motion’, J± = dQt/dt = i[H, Qt], where 
Qt is the operator of electron charge in the top layer. 
The density matrix is found from the von Neumann 
equation 





( 4 ) 


being interested in the linear response, we commute only 
the zeroth-order density matrix p^^^ with Hint- Consid¬ 
ering the harmonic time dependence, Eq. o is readily 
solved in the diagonal basis of Hq leading to a closed- 
form expression for the frequency- dependent non-local 
tunnel conductance 


Gq o; — 72 


.fsL £S 

Jp Jp' 


■'I' 


p,p'^p+q 

,s,s 


UJ ■ 


iS — [gp — 


( 5 ) 


Here, su> are the overlap integrals between the SP field 
profile and the z-components of the electron wave func¬ 
tions, sii> = (iz4'*(z)s(z)4'//(z), Upp, are the over¬ 

lap factors between the chiral envelope wave functions 
in graphene, Upp, = [1 + ss'cos0pp']/2, and are the 
occupation numbers of the eigen states of given by the 


Fermi functions shifted by the applied bias eV in the 
energy scale. 

At this point, it is noteworthy to mention the similar¬ 
ity of Eq. @ and the expressions for the graphene po¬ 
larizability 1^. The latter diverge on the ’Dirac cone’ as 
[uj^ — and a similar singularity appears in the 

tunnel conductance ([5]), except the frequency w should be 
replaced with uj — ■\/A'^ + 472^. This singularity is simply 
elucidated after extracting the real part of Eq. dSD and 
passing to the elliptic coordinates, which leads us to the 
following expression 


ReG 


J_q,a; — 


(7^s±72I 


2TTh 


{qvoY - {eVi 2 - 


( 6 ) 


where we have introduced the effective energy shift be¬ 
tween the Dirac points eVi 2 = VA'^ + 472^, and a dimen¬ 
sionless non-singular function I appearing as a result of 
the Fermi distribution integration, 


I = sinh 


eV — a;\ 

) 


f _ - ^dt _ 

J cosh -I-cosh (^4^^) 


The real part of the tunnel conductance is negative 
provided w < eV, which implies that the tunneling elec¬ 
trons rather loose their energy by emitting a plasmon 
than get the energy by plasmon absorption. At the same 
time, the tunnel transitions from top layer to the bot¬ 
tom one accompanied by SP absorption are prohibited 
by the energy conservation law as far as the SP veloc¬ 
ity exceeds the Fermi velocity. Still, just like in most 
tunnel phenomena, the negative conductivity due to the 
tunneling is proportional to a small exponent {k 

is the decay constant of the electron wave function) that 
eventually enters s± and 72. To enhance the negative 
conductance, the materials with small effective (tunnel¬ 
ing) mass TO* and small band offset Ub with respect to 
graphene have to be used. Boron nitride {Ub = 1.5 eV, 
TO* = O.Stoq) is not the best candidate for the realization 
of the net SP gain, while chalcogenides of molybdenum 
M 0 S 2 (Ub = 0.29 eV, TO* = 0.43too) and tungsten WS 2 
(Ub = 0.4 eV, TO* = 0.28too) [i^ demonstrate signifi¬ 
cantly stronger tunneling. 

The smallness of tunnel current is to a considerable 
extent compensated by a singularity in plasmon gain, 
which occurs provided eVi 2 = uj + qvo = uj{l -I- vq/s). 
The structure of the p-integral in Eq. shows that 
the singularity comes from the tunneling between elec¬ 
tron states with collinear momenta in the neighboring 
layers. Once the carrier dispersion is linear, these states 
have the same velocity and thus interact for an infinitely 
long time [3l[. The carrier scattering naturally destroys 
this coherence making the singular tunnel conductance 
finite. To account for this effect quantitatively, we add 
the imaginary self-energy corrections ImE(p,if) to the 
quasi-particle energies in These corrections appear 
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as a result of electron-phonon scattering, the latter also 
leading to the finite Drude absorption. For simplicity, 
we take ImE(p,iil) to be energy- and momentum in¬ 
dependent and equal to its value at the Fermi-surface 
ImS(pF,eF) = V, 


epTD^ 


( 8 ) 


where Z? ~ 30 eV is the deformation potential of 
graphene [i^, p = 7.6 x 10“^ kg/m^ is its mass density, 
and c ~ 0.02uo is the sound velocity 0- 

To judge on the possibility of the net plasmon gain, one 
has also to evaluate the frequency-dependent non-local 
in-plane conductivity, (Tq,ui- In the absence of electron 
collisions, one can readily apply the Kubo linear response 
theory and obtain the well-known result including inter- 
and intraband contributions (0 . The real part of in¬ 
terband optical conductivity {q = 0) is universal at high 
frequencies, Recinter = e^/4/i. However, at finite wave 
vector - which is the case of plasmons - the frequency 
dependence of conductivity becomes more complicated, 


4ft v^a;2 - q'^i 


Recrinter(q,w) ——^==== [fo i-uj/2) - fo (a;/2)]. 

( 9 ) 

The neglect of the spatial dispersion results in an under¬ 
estimate of Reuinter and, hence, plasmon damping. This 
is especially crucial for the acoustic modes which velocity 
just slightly exceeds the Fermi velocity and the respec¬ 
tive interband conductivity is in a dangerous vicinity of 
the singularity at w = qvQ. 

Considering the intraband (Drude) conductivity, a spe¬ 
cial care should be taken to account correctly both for the 
spatial dispersion and finite carrier relaxation time [l6| . 
Apparently, the effects of carrier collisions cannot be in¬ 
cluded by a simple replacement a; —>■ a; + w in the Kubo- 
like expression for the collisionless conductivity (i^ . 
Even in the quasi-classical formalism of Boltzmann equa¬ 
tion (which is justified in the case of interest, q qp, 
uj <C ep), the collision integral in the r-approximation 
violates the particle conservation. To avoid those diffi¬ 
culties, we have evaluated the conductivity by solving the 
kinetic equation with the particle-conserving Bhatnagar- 
Gross-Krook collision integral in the right-hand side. 
The respective intraband dynamic conductivity reads 


^intra(q5 — 


ie^sp 


^ ^ UJ ^ UJ ■' 


^_^ u qvQ J / qvQ \ ’ 

27r UJ UJ ^ 


( 10 ) 


where Jn{a,b) = dd cos^ 9 /[I — a cos 6 + ib], ip = 

Tln(l-|-e^^/^), and p is the collision frequency limited by 
the electron-phonon scattering 0 given essentially by 
Eq. ([H). At zero wave vector, Ji = 0, and we restore the 
Boltzmann conductivity crintra(0,w) = ie^ipjTr{uj + iv), 
while for nonzero q the real pat of intraband conductiv¬ 
ity typically appears to be larger than its zero-g value. 


This difference between no-local and local conductivities 
agrees with the experimentally observed distinction be¬ 
tween the plasmon lifetime and the Boltzmann relaxation 
time 0. 



FIG. 2. Net ’’effective conductivity” as it appears in the 
SP dispersion law, Re[crq,;j -|- 2Gq^uijq'^] (in the units of e^/ft) 
vs. frequency and wave vector at two different temperatures: 
T = 77 K (top) and T = 300 K (bottom). The barrier layer 
is 2.5 nm WS 2 . The singularities of the interlayer tunnel 
conductance, hw = eVi 2 ± gwo are shown with blue lines, the 
dispersion of acoustic SP is shown with red line. In the region 
filled with green, the net effective conductivity is negative. 

With these prerequisites, we are able to compare the 
plasmon gain due to the tunneling and loss due to 
the Drude and interband absorption. Namely, if in 
some frequency range the ’effective conductivity’ (Jes = 
Re[(Tinter + Sintra — 2G/g^] < 0, then the plasmon gain ex¬ 
ceeds loss. The self-excitation of plasmons with frequen¬ 
cies satisfying the net gain condition is possible. Figure^] 
shows the wave vector- and frequency dependence of the 
effective conductivity at the room and nitrogen tempera¬ 
tures with the region under the green plane correspond¬ 
ing to the negative values. The two sets of singularities 
are clearly present in the conductivity spectra: the ab¬ 
sorption singularity at a; = qvo comes from the in-plane 
conductivity, and the gain singularity at w = eVi 2 ± qvQ 
comes from the enhanced tunneling between the collinear 
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FIG. 3. Comparison of the real part of in-plane graphene con¬ 
ductivity Re[crinter + cTintra] governing the plasmon loss (nor¬ 
malized to jh, dashed lines), and the real part of the ef¬ 
fective tunnel conductance —2ReG/g^ governing the plasmon 
gain (solid lines). The barrier layer is 2.5 nm WS 2 . In the fre¬ 
quency ranges where Re[(Tinter + crintra — 2G/q^\ < 0, the plas- 
mons are amplified instead of being damped. These ranges 
are marked with color bars near the frequency axis. 


states. As the wave vector approaches zero, the effective 
conductance grows indefinitely due to the presence of 
in the denominator, this will not, however, lead to the 
infinite plasmon gain as the acoustic SP spectra develop 
a gap Aw ~ n at small wave vectors On the other 


hand, there is a wide range of frequencies w > O where 
the dispersion of acoustic SPs in not strongly affected by 
tunneling, but the real part of net effective conductivity 
is negative, which justifies the possibility of gain. 

The spectra of in-plane conductivity Re[(Tinter + Sintra 
and effective tunnel conductance — 2Re[G/g^] are com¬ 
pared in Fig. [3] for the wave vectors satisfying the acous¬ 
tic plasmon dispersion w = sq. At room temperature 
(upper panel) the self-excitation is possible in a narrow 
(^ 15 — 20 meV) vicinity of the tunneling resonance fre¬ 
quency, while at higher frequencies the strong interband 
absorption surpasses the gain. At lower temperature of 
77 K (lower panel) the range of frequencies correspond¬ 


ing to the net SP gain is significantly broader, while the 
peak tunnel conductance is higher. The broadening of 
the gain region is due to efficient low-temperature Pauli 
blocking of the interband transitions with energy below 
the double Fermi energy. The sharpening of the reso¬ 
nant peak is attributed to the phonon freezing-out at 
low temperatures, which increases the quasiparticle life¬ 
time and resonant gain. It is also noteworthy that the 
peak in the negative tunnel conductance occurring at 
Wres = -I- uo/s) ~ eRi 2/2 falls into the ’trans¬ 

parency window’ of graphene where the Drude absorp¬ 
tion is low (wres ^ 1^) but the interband absorption is 
still blocked (wres < 

The efficiency of SP excitation can be further improved 
in the gated double-graphene layer structures, where the 
energy shift between Dirac points A and the carrier den¬ 
sity can be controlled independently. At the fixed carrier 
density, the plasmon gain increases with reduction in A 
as the energies of the ^ = -1-1 and I = —1 states get closer 
in energy scale resulting in enhanced tunnel coupling. 

In realistic graphene-based tunnel structures, the ad¬ 
jacent layers are always slightly misaligned in real space, 
which results in the misalignment of Dirac cones in the 
A:-space. The finiteness of SP wave vector can, to some 
extent, compensate this misalignment. If the Dirac cones 
in adjacent layers are separated by a wave vector q^/, the 
tunneling resonance will be retained, though the position 
of the resonant peak will depend on the direction of plas¬ 
mon propagation q. The resonant condition in this case 
is easily shown to be 

|q + qMpu^ = (eFi2 - a;)^ (11) 

and the expressions for the real part of tunnel conduc¬ 
tance are obtained from Eq. ([S]) by a simple replace¬ 
ment g — |q + qA/l- The case of the plasmon-assisted 
resonant tunneling is thus radically different from the 
photon-assisted tunneling (^ . where even a slight mis¬ 
alignment breaks the resonance due to the negligible pho¬ 
ton momentum. The robustness of SP emission against 
slight misalignment correlates with the broadness of gain 
(and emission) spectra extending above the resonant fre¬ 
quency. At the same time, the spectra of phonons emit¬ 
ted upon tunneling are predicted to be Lorentzian their 
width being limited by the carrier relaxation rate. For 
these reason, the inelastic tunneling accompanied by the 
emission of surface plasmons rather than photons may 
largely contribute to the spontaneous terahertz emis¬ 
sion from the double graphene layer structures observed 
in [^ . The mechanism of electromagnetic radiation in 
case of plasmon excitation could be either the dipole radi¬ 
ation of the whole heterostructure, or the radiative decay 
of plasmons scattered by the edges. 

In conclusion, we have theoretically studied the 
frequency-dependent non-local tunnel conductivity in the 
biased graphene layers. In a wide range of frequencies, 
the real part of the tunnel conductivity is negative due 
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to the dominance of inelastic tunneling with surface plas- 
mon emission over the absorption. Moreover, at frequen¬ 
cies satisfying the condition huj = eVi 2 ± guo, where V 12 
is the interlayer potential difference and q is the plas- 
mon wave vector, the absolute value of negative tunnel 
conductivity is resonantly large, which is a manifesta¬ 
tion of enhanced tunneling between electron states with 
collinear momenta in the neighboring layers. A detailed 
comparison of the surface plasmon loss due to the inter¬ 
band and Drude absorption and gain due to the inelastic 
tunneling shows the possibility of the net gain in suffi¬ 
ciently thin tunnel structures. 

The work of DS was supported by the grant # 14- 
07-31315 of the Russian Foundation of Basic Research. 
The work at RIEC was supported by the Japan Soci¬ 
ety for Promotion of Science (Grant-in-Aid for Specially 
Promoted Research No. 23000008). 

APPENDIX 

Evaluation of the in-plane conductivity 



Energy, ftco (meV) 

FIG. Al. Real part of the interband in-plane conductivity 
vs. frequency calculated using exact expression (IaHi and 
its various approximations: local (q = 0), black dashed line; 
Eq. (IMI, green dashed line; Eq. (lA7l) . red dashed line. For 
the non-local results, the wave vector q satisfies the dispersion 
of acoustic plasmon q-{ui) = uo/s. Temperature T = 300 K, 
Fermi energy £f = 75 meV, s = l.luo- 


The expression for the interband part of conductivity 
is readily obtained from the Kubo theory (45| 


Winter (Qi — 



E 



ss' 

PP 

2 

£s' £S 

Jp' Jp 


e®) — e® 

p' p 

[w -b 1(5 - (e®' - e® )] 


(Al) 


Here Vp_^ p_ is the matrix element of velocity operator 
in graphene v = vqa, p± = p ± q/2, /* is the electron 
distribution functions in the valence (s = —1) and the 
conduction (s = -|-1) band, and ep is the dispersion law 
in the s-th band. Known the eigen functions of graphene 
Hamiltonian Hq = r'oO’Pi 


|sp) 


1 

72 



(A2) 


one readily finds (cp_ | |r;p+) = 

ivq sin [(0p+ + 0p_)/2]. The subsequent calculations are 
conveniently performed in the elliptic coordinates 


p = - {coshitcosr;,sinhMsinu} . 


(A3) 


In these coordinates |p±| = (( 7 / 2 )[coshu ± cosu], 

\ {cp-\vx \vp+) I'^dpxdpy = {qvo/2)"^ cosh^ usin^ vdudv. 
The real part of conductivity is readily extracted with 
Sokhotski theorem, which leads us to 


RefJinter(q: ^) — 


W 


\JIjP' — q 


2^,2 


dv sin"^ V X 


|/o 


w qvo 


■ cos V 


- fo 


| + ^cos?;j|. (A4) 


To proceed further, we note that in the domain of inter¬ 
est w > qvo one always has qvq cosv < w. As a simplest 
approximation, one can set 5 = 0 in the arguments of 
distribution function and obtain 

2 

Recri„ter(q,w) Ri ^—=^== [/q (-w/2) - /q (a;/2)]. 

(A5) 

The agreement between approximate analytical expres¬ 
sion (IA5|) and exact integral relation (IA4|) is fair as far 
as qvq is not too close to w. A better agreement can 
be obtained after the change of variable cos u = t iw 
Eq. (IA4I1 by noting that sinu = \/l — varies rapidly 
over the integration domain t S [— 1 ; 1 ], while the re¬ 
mainder g{t) = /o[—(w — qt)/2] — /o[(w -I- qt)/2\ varies 
slowly. Thus, one can make the following approximation 

1 11 

J g(t)'\/r^7^dt Ri - J \/l — t'^dt j g{t)dt. (A6) 

-1 -1 -1 

The following approximation for the conductivity holds 


Re(Tinter(q; ^) 


T uj 

_^ 

4/1 qvo _ q 2 y 2 

cosh ^ -I- cosh 
cosh ^ -I- cosh 


(A7) 


Figure lATl shows the calculated interband conductiv¬ 
ities using the exact jAil and approximate expressions 
(IA5I) and (IAtI) as well as the local {q = 0) approxima¬ 
tion. Clearly, the neglect of spatial dispersion in the case 
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of acoustic SPs with velocity slightly exceeding the Fermi 
velocity results in an underestimation of the damping. 

We now pass to the in-plane conductivity associated 
with the intraband transitions, Re(Tintra(qj w). We shall 
restrict ourselves to the classical description of the in- 
traband electron motion which is justified at frequencies 
hiu £f, q ^ <lF- The effect of plasmon tunneling gain 
occurs in this frequency domain, otherwise, the inter¬ 
band damping of SPs takes place [see Eq. ig. To account 
for the carrier relaxation in the non-local {q ^ 0) case we 
adopt the formalism of kinetic equation with the particle- 
conserving Bhatnagar-Gross-Krook collision integral T71 
in the right-hand side. 


- *W(5/q(p) -h iqv(5/q(p) -H ieqytpq-^ = 
xf t \ dfo 


— V 


(AS) 


Here <5/q(p) is the sought-for field-dependent correc¬ 
tion to the equilibrium electron distribution function /o, 
i5nq is the respective correction to the electron density, 
V = uqp/p is the quasi-particle velocity, and v is the elec¬ 
tron collision frequency which is assumed to be energy- 
independent. The current density, associated with the 
distribution function (5/q(p) reads: 


l^jq = 


-egY. 


P 


^ d/o zegvyjq - w^Sn^ 

de u! + iv — qv 


(A9) 


Recalling the relation between small-signal variations of 
density and current, ujSuq = gdjq, and evaluating the 
integrals in Eq. (Ei, we find the conductivity given by 

Eq. uni). 


Spectra of plasmons coupled to the double graphene 
layer heterostructures 


The plasmon spectra are obtained by a self-consistent 
solution of the Poisson’s equation 


- q^5^pc^{z) 
47r 


dz"^ 


-[i5pt,q(5(2: - d/2) + Spb,qS{z + d/2)], (AlO) 


the continuity equations 


- iwSpq -I- qdjq = 0, (All) 

and the linear-response relation between current density 
and electric field, (5jq = Uq^^dEq. Here q is the two- 
dimensional plasmon wave vector, d is the distance be¬ 
tween layers, n is the background dielectric permittivity, 
5pt,q and dpt,,q are the small-signal variations of charge 
density in the top and bottom layers, respectively. In 


the absence of built-in voltage, due to the electron-hole 
symmetry, the charge densities in the layers are equal in 
modulus an opposite in sign, moreover, the layer conduc¬ 
tivities are equal. This allows us to seek for the solutions 
of Eq. (lAlOl) being symmetric and anti-symmetric with 
respect to the electric potential. A straightforward cal¬ 
culation brings us to the following dispersions [35 1 


— ( 1-f cothy ) -f guq^ =0 (A12) 


for the antisymmetric (acoustic) mode, and 


- ^ + tanhy ) + gcTq,^ = 0 (A13) 


for the symmetric (optical mode). Being interested in 
the long-wavelength limit, qd/2 <C 1, we perform the 
expansions 1 -I- cothqd/2 k, 2/qd, 1 -|- tanhqd/2 Ri 1. 
In the same limit, the conductivity is essentially classi¬ 
cal, moreover, the interband transitions do not affect the 
low-energy part of the spectra. With these assumptions, 
we use the following (collisionless) approximation for the 
conductivity 




e sf ^ 


h 2 tt q'^VQ 






- 1 


(A14) 


which follows readily from Eq. (US after setting f = 0. 
Equation (IA12|) admits an analytical solution, which is a 
sound-like dispersion 


OJ- 


Vo 


1 -I- AacqFd 
v/1 + SacqFd 


(A15) 


Here, we have introduced the Fermi wave vector qF = 
£f/vq. The velocity of the acoustic mode always exceeds 
the Fermi velocity thus preventing the Landau damping. 
The dispersion equation for the optical mode io+{q) is 
cubic, however, in the long-wavelength limit the spatial 
dispersion of conductivity can be neglected as the phase 
velocity of this mode significantly exceeds the Fermi ve¬ 
locity. The approximate relation for Ct)+(q) has the fol¬ 
lowing form 


UJ+ pz vo^JAacqqF■ 


(A16) 


To account for the damping of the propagating plas¬ 
mon one can express the wave vector q from Eqs. (IA12I) 
and (IA13I) and consider the real part of conductivity 
along with the imaginary one. The results of such cal¬ 
culations are shown in Fig. IA2I along with the frequency 
dependencies of the net conductivity. The acoustic mode 
is typically damped more heavily than optical one the 
two effects being responsible for that fact. First of all, 
the velocity of optical mode is higher than that of acous¬ 
tic one, hence, at a fixed scattering time, the free path of 
optical mode is higher. The second reason is the increase 


























FIG. A2. Top panel: real part of the non-local net in-plane 
conductivity (Drude -|- interband) vs. frequency calculated 
for the wave vectors satisfying the dispersion of the acoustic 
mode q(Uj) (orange line) and the optical mode q+{uj). Tem¬ 
perature T = 77 K, Fermi energy ef ~ 150 meV, s = 1.2uo. 
Bottom panel: calculated imaginary parts of the propagation 
constant for the acoustic and optical modes calculated using 
the local and non-local expressions for the conductivity 


mon [Hq in Eq. ((2)] constitutes the blocks describing 
isolated graphene layers Hg± and the block describing 
tunnel hopping T. This Hamiltonian acts on a four- 
component wave function 

if = (A19) 

which components are the probability amplitudes of find¬ 
ing an electron on a definite layer and on the definite lat¬ 
tice cites, A or B. Such description of electron states is 
common for graphene bilayer; moreover, the elements of 
tunneling matrix have been evaluated for different twist 
angles between single layers constituting a bilayer. In 
more comprehensive theories, the T-matrix is affected by 
the band structure of dielectric layer. Here, for the sake 
of analytical traceability, we choose the tunneling matrix 
in its simplest form which is applicable to the AA-stacked 
perfectly aligned graphene bilayer, 'T = HI, where H can 
be interpreted as the tunnel hopping frequency. 

To estimate its value, we switch for a while from the 
tight binding to the continuum description of electron 
states in the z-direction. We model each graphene layer 
with a delta-well U{z) = Ad(z — Zt^b), where A is the po¬ 
tential strength chosen to provide a correct value of elec¬ 
tron work function Ub from graphene to the background 
dielectric, 


A = 



(A20) 


where m* is the effective electron mass in the dielec¬ 
tric. The effective Schrodinger equation in the presence 
of voltage bias A/e between graphene layers takes on the 
following form 


in the real part of conductivity at high spatial dispersion 
discussed above. The use of ’local’ expressions for con¬ 
ductivity results in a factor of two error in the free path 
of acoustic mode, which is also illustrated in Fig. ( IA2I) . 

In the following calculations we shall also require the 
spatial dependence of the plasmon potential in the acous¬ 
tic mode, which can be obtained from (lAIOl) . It is con¬ 
venient to present it as 


Sipq(z) = S(pos(z), (A17) 

where <po is the electric potential on the top layer. The 
shape function has the following form 

g-fc(z+d/2)^ z < -d/2, 

sinh {qz/2) 


siz)= { - 


|z| < d/2, 


(A18) 


sinh {qd/2) 

_^-k(z-dl2)^ z > d/2. 

Estimate of the effective tight-binding parameters 


The tight-binding Hamiltonian of the tunnel-coupled 
graphene layers in the absence of the propagating plas- 


+ C/j’(z)^(z) = A4'(z), (A21) 

where C/f is the potential energy created by the applied 
field 


r 1, z < -d/2, 

Up (■^) = y 'I 2-z/c^, l-^l < d/2, (A22) 

[ - 1, z > d/2. 

The solutions of effective Schrodinger equation represent 
decaying exponents at |z| > d/2, and a linear combina¬ 
tion of Airy functions in the middle region |z| < d/2 

^niz) = CAi {—z/a + e) + DBi {—z/a + e), (A23) 

where £ = 2m*\E\a^/fi^ is the dimensionless energy and 
a = {h‘^d/2m*Ay/^ is the effecive length in the electric 
field. A straightforward matching of the wave functions 
at the graphene layers yields the dispersion equation 
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det 


/ g-kid/2 

{2kb — ki) 

0 

V 0 


—Ai {d/2a + e) 

— iAi' ((i/2a + e) 
—Ai {—d/2a + e) 
— iAi' {—d/2a + e) 


where kb = y^2m*Ub/fi^ is the decay constant of the 
bound state wave function in a single delta-well, fci = 
y/2m*{E + Al2)lh?, k 2 = y^2m*{E — Aj2)jh? Equa¬ 
tion (IA24|) yields two energy levels Ei {I = ±1) which 
cannot be expressed in a closed form via the parameters 
of the well. However, the dependence of Ei on the energy 
separation between layers A is accurately described by 


EiiA) = -C/,+1 ^(E+i,a=o - A_i,a=o)" + (A25) 


The energy spectru m (IA25I) is typical for the tunnel cou¬ 
pled quantum wells ^ ; the same functional dependence 
of energy levels on A is naturally obtained by diagonal¬ 
izing the block Hamiltonian ([2]), 


EiXA) = -Ub + (A26) 

This allows us to estimate the tunnel coupling H as half 
the energy splitting of states in double graphene layer 
well in the absence of applied bias 

= 2 [^-1-1,A=o “ A_i^a=o] • (A27) 

The Z-index governs the ^-localization of electron in a 
biased double quantum well. At large bias A H, the 
delta-wells interact weakly, thus Z = -|-1 corresponds to 
the state localized almost completely in the top layer and 
Z = — 1 corresponds to the electron in the bottom layer. 
At small bias A Ri H the state Z = -1-1 is anitisymmetric 
and Z = — 1 is symmetric. 


Electron-plasmon interaction and solution of the von 
Neumann equation 

The presence of plasmon propagating along the double 
graphene layer results in an additional potential energy 
of electron 


i7i„t(r,t)=e^(z)e*(«"--‘), (A28) 

where we assume the direction of plasmon propagation 
to be along the cc-axis, and the dependence of potential 
on the z-coordinate is given by Eqs. (IA17I) and (IA18I) . 
The additional terms in Hamiltonian due to the vector- 
potential are negligible as far as the speed of light sub¬ 
stantially exceeds the plasmon velocity. 

With our choice of the tight-binding basis functions 
(IA19|) as those localized on a definite layer and on a def¬ 
inite lattice cite, we shall require 16 matrix elements of 


—Hi {d/2a E e) 

— iBi' {d/2a E e) 
—Hi (—cZ/2a -|- e) 
— iBi' {—d/2a E e) 


0 \ 

0 

g-/c2d/2 


= 0, (A24) 


the potential energy (IA28I) connecting those basis states. 
However, it is more convenient to work out the matrix 
elements of (IA28I) connecting the eigen states of Hamilto¬ 
nian The good quantum numbers of these states are 
the in-plane momentum p, the band index s = ±1 (-1-1 
for the conduction band and —1 for the valence band) 
and the Z - index discussed above. The respective matrix 
elements are 


(p, s, Z|i7int Ip's'Z') = 

^p,p'-q«PP'e7’o / (A29) 


Having obtained the matrix elements of electron- 
plasmon interaction, we pass to the solution of the von 
Neumann equation ®. Being interested in the linear re¬ 
sponse of electron system to the plasmon field, we decom¬ 
pose the density matrix as p = E where p^^'> is 
linear in the electron-plasmon interaction i7int. The cal¬ 
culations are conveniently performed in the basis of the 
eigenstates of Hq including the effects of tunneling non- 
perturbatively. In this basis, [Hq, p^^^ad = i^a~ 
and thus one readily writes down the solution for the 
density matrix 


Pap 


.^int) P 


(0) 


uj E id — {Sa — ep) 


(A30) 


The first-order correction is now expressed through the 
density matrix in the absence of plasmon field 
A particular choice of p*^°^ requires the solution of ki¬ 
netic equation in the voltage-biased tunnel-coupled lay¬ 
ers, however, in several limiting cases the situation is 
greatly simplified 28[. If the tunneling rate H is slower 
than the electron energy relaxation rate (e.g., due 
to phonons and carrier-carrier scattering), the quasi¬ 
equilibrium distribution function is established in each 
individual layer. In this situation, p^°^ is diagonal in 
the basis formed by the wave functions localized on top 
and bottom layers its elements being the respective Eermi 
distribution functions. In the other limiting case, when 
tunneling is stronger than scattering (H ^ p^), the elec¬ 
tron is ’collectivized’ by the two layers, and the density 
matrix p^^'> is diagonal in the basis of Ho-eigenstates. Eor 
the parameters used in our calculations. Ml r: 10 meV 
exceeds the relaxation rate Zip r 1 meV, and the latter 
limiting case is justified. Setting p^^^j = faSap, where / 
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is the Fermi distribution function, we find 
\p's'l') = 


£s'l' _ £sl 

OJ + 10 — [Ep — Ep, ) 
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